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1 Introduction 



This paper is concerned with a numerical method for the Cauchy problem 



u{x, 0) = /(x), X G 



for exponents m > 1 and space dimension > 1. We present a numerical 
method for this equation. We also prove existence and uniqueness of solution 
to the method, via a maximum principle. Moreover, convergence to the 
theoretical solution is also proven. We extend later the results for equations 
with (p{u) instead of (in the following u"^), for some monotone (f 

with good regularity conditions. The general theory of existence, uniqueness 



and regularity of solutions for the equation (2.1) has been studied by A. de 
Pablo, F. Quios, A. Rodriguez and J.L. Vazquez in [6]. They also study a 
more general case with (— A)"^/^ for a G (0,2) in [7]. Even more, in pj they 
study the case where the diffusion is logarithmic, that is, ip{u) = \og{u + 1) 
as de natural limit as m — of </?(«) = u"^. 

We recall that the nonlocal operator (— A)^/^ is well defined via Fourier 
transform for any function / in the Schwartz class as the operator such that 

or via Riesz potential for a more general class of functions, 
(-A)'/VW = C.RV, / /f 

Jrn f ~ y\ 



N+l . 



where Cn = ^ ^{^^) is a normalization constant. For an equivalence 



of both formulations see for example [TTj . 

Previous works in numerical analysis for nonlocal equations of this type 
are done by S. Cifani, E. R. Jakobsen, and Karlsen in [2], [3], [1]. In particular 
they formulate some convergent numerical methods for entropy and viscosity 
solutions. One of the main differences of the present work is that we don not 
deal directly with the integral formulation of the fractional laplacian, instead 
of this, we pass through the Caffarelli-Sylvestre extension ([Ij) with implies 
solving a problem with only local operators in one more space dimension. 
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2 Local formulation of the non-local problem 
2.1 The problem in 

Our aim is to find numerical approximations for tlie solutions of the next 
porous medium equation with fractional diffusion, 



(2.1) 



du 
di 

u{x,0) = f{x) 



(a;,t) + (-A)i/2^™(x,t) = xeR^,t>0, 

N 



X E 



with m > 1 and / G L^(]R^) nonnegative. The general theory for existence, 
uniqueness and regularity of the solution of problem (2.1) can be found in 



[H]. In particular, they state that problem (2.1) is equivalent to 



( Aw{x,y,t) =0, 



(2.2) 



X G M^, y>0, t>0, 



-^{x,0,t) = —{x,0,t), xgM^, y = 0, t >0, 
L w{x,0,0) = /™(x), X G M^, 

in the sense of trace and harmonic extension operators, that is, 

u{x, t) = Tr(w^/"'(x, y, t)), w{x, y, t) = E{v^{x, t)). 



In (2.2), A denotes the + 1 dimensional laplacian operator, 

A- + ^ 
■p^ dx} dy'^ 

2.2 The problem in the bounded domain 



In order to develop a numerical solution to Problem (2.2), we are going to 



perform a monotone approximation of the solutions for the problem posed a 
bounded domain. 

We consider X;, F, T G M positive for / = 1, A^. We define the bounded 
domain Vt = (— Xi,Xi) x ... x {—Xn,Xn) x (0,y), and T = d^l. For con- 
venience we also divide the boundary in two parts: = [— Xi,Xi] x ... x 
[—Xn,Xn] X and Th = 9f2\Frf. With this notation, we formulate the 
problem in the bounded domain as 
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Aw{x,y,t) = 0, 



(2.3) 



{x,0,t) 



dw 



dt dy 
w{x,y,t) = 0, 



{x,y)eQ, te{0,T], 

x,o,t), {x,y)era, te{o,T], 

{x,y) G r^, 
{x,y) G Th, 



where we have imposed homogeneous boundary conditions for Th- 
in the following we will consider the problem with = 1 in order to 
simplify de notation but all the arguments are also valid for > 1 without 
any extra effort. 



Y 



+ 

I 
I 



n 



-X 



II 

X 



3 Discrete formulation 



In order to solve problem (2.3) fot t G [0,T], first we perform a time and 
space discretization. For time discretization we choose the number of steps 
J, and then 

< J At < T, j = 0, J where At = T/J and tj = jAt. 

We also need to discretaze the space domain fl = [—X,X] x [0,F]. Let 
/, K be the number of steps on each space direction, 

< iAx <2X, i = 0, / where Ax = 2X/I and Xi = iAx - X, 

< kAy <Y, k = 0,...,K where Ay = Y/K and yk = kAy, 
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Vk 

2/2 

yi 

2/0- 



Xq Xi X2 ••• Xi 



Xi 



We use the notation 
(3.1) 



w{xi,yk,tj) = (wj) 



for the value of the solution w to Problem (2.3) in the points of the mesh, 
and 

(3.2) w{x,,yk,t,) ^ {Wj)t 

for the solution of the numerical method. 



3.1 Numerical Method 

Assume that Ay = Ax. On each time step j = 1, J we have to solve the 

following linear system of equations 

(3.3) 

( W)JV, + TOt. + TO?«+iMl^i(M.O. 0<,<I.O<k<K, 



Ax^ 



At 
Ax 



(w,)1 = 0, 



where the solution of 

r (^o)ti + (H^o)ti + (W^o)^^ + {Wo)t' - mo)' 

(Wo)', = r (x.), 

I iWo)1 = 0, 
is used to start the numerical method. 



, if < 2 < J, 
otherwise, 

0, <i < 1,0 < k < K, 
if < z < /, 
otherwise. 
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3.2 Local truncation error 

We will say that the local truncation error {Tj)f is the error that comes from 



plugging the solution w to Problem (2.3) into the numerical method (3.3). 
Let us also write 



(3.4) A = max|(rj 



Theorem 3.1. Let w be the solution to Problem (2.3). Assume that there 
exist two constants Ci,C2 > such that 

CiAx <At< C2AX. 

Then, 

(3.5) A = 0(At{Ax + At)). 

Proof. Of course the local truncation error in the boundary nodes situated on 
the part Th of the boundary is zero since we have imposed that the solution 



is zero in Th and equal to f"^{x) in Yd as in Problem (2.3) 



If < z < / and < k < K (the interior nodes), then 
1 



= Aw{xi,yk,tj^i) + O(Ax^) = 0{At{Ax + At)). 
If < z < / and k = 0, then 

dw dw^^^ 
= At[—{xi,0,t,_i) + O{Ax)] -At[^^(x„0,Vi) + O(At)] 

= 0{AtAx) + 0{At^) = 0{At{At + Ax)). 

□ 

3.3 Existence and uniqueness of the numerical solution 

We define now the quantity needed for the maximum principle theorem, 

bmax = max{/'"(a;)}. 

X 

In the following we will denote by (f{x) = x^. If m > 1 then '^'{x) = mx™~^ 
is a locally bounded function for a; > 0. 
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Theorem 3.2. [Discrete maximun principle] Let iWj)^ he the solution to 
Problem (3.3) with m>l. and assume At < C{m, f)Ax where 

(3.6) C{mJ) = [m{b^a.Y"'-'Y'- 
Then, for every i,k,j we have 

(3.7) < (W,)^ < b^a.. 

Note 1. It is interesting to remark that C{l,f) = 1, which means that we 
recover the expected restriction At < Ax for the linear case. In the literature, 
this kind of condition use to be called CFL condition. 

Proof. On each time step we have a discrete harmonic extension problem and 
so, it is sufficient to prove the maximum principle in the boundary nodes and 
therefore the interior nodes are automatically smaller than them. We will do 
the proof by induction on each time step. It is trivial that 

< (Wo)^ < b^a.. 

Now we assume that 

< iW,^i)1 < b^a.. 

Then, 

and changing variables to (f/j)f = [(W^j)i^]"'^^™', we obtain that 

At 



(3.8) (f/,)," = ^im-i)]r - m,-im + (t^.-i)°- 

Now, using the Mean Value Theorem 



for some ^ G [{Uj-i)}, (f/j-i)^]. Then we can rewrite (3.8) as 



(3.9) 



1 - ^'(0 



At 



l(f/.)°l < v'{0^\iu,-i)l\ + 



Ax 

At this mome nt , thanks to our induction hypothesis and the value of the 

constant (3.6), it follows that — — ''^^"^'^^ — ^ 1 therefore 
J As Ax 

At 



1 - ^'(0 
1 - ^'(0 



Ax_ 
At 
Ax 



l(f/.-i)° 
{b. 



l/m 



{br. 



The same argument holds for (f/j)° > 0. 



□ 
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Corolary 3.3. If At < C{m, f)Ax, then Problem {3.3) has a unique solu- 
tion. 

Proof. We are only going to proof uniqueness, and since we are working with 
a linear system of equations, the existence is equ ivalent to the uniqueness. 
Let {Wj)'l and (Vj)f be two solutions of (3.3) and let us define 



Then (lo)f satisfies (3.3) with / = and so, by the discrete maximum 
principle, 

< (lo)- < 0, 

for all z, k. Proceeding by induction we get that for all z, j, k we have 

(y^t = 0. 



□ 



3.4 Convergence of the numerical solution 

Since we are originally interested in the solution at the boundary {Uj)^ = 
[{Wj)^Y^'^ we have two options in order to define the define de error of the 
numerical method. The first option, 

(3.10) if,)^ = wix,,y,,t,)-{W,)l F,=max|(/,)f|, 

t,k 

and the second one 



(3.11) (ej)f = u{xi, Vk, tj) - {Ujfi, Ej = max |(ej 

i,k 



But anyway, if we are able to control (3.11) we have also controlled (3.10) 
because ip'{x) is locally bounded and 



= KOf - TOf = [Mr - [(^.) 
= [(«.H-(f^)-]¥''(0 



for some ^ G [{ujYi, {Uj)f]. This implies that |(/j)f| < C(m, /)|(ej)f | and 
therefore Fj < C{m, f)Ej. 
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Theorem 3.4. Let w be the solution to Problem (2.3) and {Wj)\ be the 



solution to system (3.3) with m > 1. Assume that there exists two constants 



C{m, f), D > such that 



Then 



for j = 1,...,J. 



DAx <At< C{m, f)Ax. 



Fj = 0{Ax + At), 



Proof. As in the local truncation error, the election for the boundary condi- 
tions in Th in the numerical method give us error zero there. 

Lets us denote {Eb)] and {Ei)j the maximum errors in the boundary 
nodes and in the interior nodes at time j'At, that is 

{Eb)] = max |(ej)°|. 

Since we have chosen a second order approximation for the laplacian, it is 
clear that 

Ej = max{{EB)j,0{Ax^)} < {Eb)j + 0{Ax^). 
If < i < /, we have the equations 



Rewriting the above equations en terms of u and {Uj)f we get 



At 
Ax 



Ax 



Subtracting them, we obtain 
(3.12) 



Ax 



-K 



By mean value theorem, relation (3.12) turns into 
At 



Ax 



(K-i);-(f^.-i).>'(6)-((«,-i)°-(t/,-i)i')v''(eo) 



+ e,- 
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for some e [(wj-i)-, (^j-i)-] and e (^j-i)J]- Then 



At 
Ax 



+ - (r,-i)°, 



that is, 

= |^^'(6)(e,-i); + [l - l^^'(eo)] (e,-i)° - (r,-0°. 
By assumption, all the coefficients that comes with (ej_i)f are positive, then 



1 - I^^Uo)' 



lfe-.)?l+A 



(3.13) 



1 - 



^,-1 + A. 



So the key point is to be able to control the difference between </?'(^i) and 
(/7'(^o) and this is not difficult at all. Since there exists a constant (assuming 
enough regularity of the solution u) K >Q such that 



■)]-{u^f,\<KI^x and \{U ^)\ - {U < K I^x . 



then (This formula is valid only for m natural, and in this case <^'(^i) 

TtA{u,)\r-'m)w^-' and = T.iU[{u,)vr-%u,)^^-') 

m 

1=1 

m 



1=1 

m 



1=1 



= {Kj^i)1 + RAx. 

But it is not difficult to prove that in fact |v?'(Ci) — < RAx where 

> is a constant depending only on m, and hmax- The proof is only 
based in the idea that, the function 



rj.m _ym 

x-y 
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is C\{0, oo) x^xoo)), and ^.'(^0) = (^,)°) and ^'(^i) = 9{{u,)}, {Uj)}). 

Ej^i + A 



Then from (3.13) we obtain 

At 
Ax 
At 



¥^'(^0) + RAx 



1 - 



Ax 



RAxEj.i + Ej^i + A 



< (1 + i?At)Ej_i + A. 
Remember also that we have Ej < max |(ej)^| + O(Ax^) and A = O(Ax^). 

0<i<I 

Then we have the next recurrence equation for the error 
(3.14) Ej < (1 + RAt)Ej^i + A. 

We should also remember that At = T/J where T was the final time and 
J the number of elements in the time discretization. Then, we can rewrite 
( |3l4l ) as 

1 

J' 

for some constant C > 0. 

Of course it is enough to bound Ej to ensure the convergence of the 
method. Since A = O(At^) lets say that there exists a constant L > such 
that A < LAt^. Then, from (3.13) we obtain 

Ej < il + Cj)Ej^, + LAt^ 



(3.15) 



Ej < (l + C-)E,-_i + A, 



< (1 + C 

< (1 + C 



1^ 
J' 

J' 



< 



<{i + c-y 



Ej^i + LAt^ 



[1 + Cj)Ej_2 + LAe + LAe 



Ej_2 + 2LAt^ 
1 



Eq + LJAr 



But (1 + C^y < e^, J At = T and Eo < DAt^ for some D >, so 



that is 



Ej < DAt^ + LTAt , 
Ej = 0{At) = 0{Ax + At). 



□ 
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3.5 Properties of the scheme 



We present now some properties that can be extracted directly from the 
numerical scheme ( |3.3[ ). They are the analougus of some of the energy stim- 
ulates presented in [6] and [7j. The first one is a direct consequence of the 
maximum principle. 

Corolary 3.5. [Comparison principle] Let f,g & C^([— X, X]) nonnegative 
such that f{x) > g{x) Wx G [— X, X] and let also (Wj)^ and {Zj)^ be the 
correspondent solutions of Problem (3.3) with m > 1. If At < C{m, f)Ax 
then, 



Proof. Let {Hj)'^ be the solution of Problem (3.3 ) with the nonnegative initial 
data h 



Then by the Discrete Maximum Principle 



3.2 



(Hj)^ > 0. We recall that at time j = the scheme is linear and so (Wq 
(Zo)f = (-f^o)f > 0. Proceeding by induction, assume that {Wj)^ - {Zj)^ > 
or equivalently {Uj)f — (V^)f > 0. The next two equations holds 

At 
' Ax 

^o_ 
- Ax 

subtracting them and using the mean value theorem 

At 

Ax 
At 



+1J 



At 



Ax 



1 



Ax 



And thanks to our CFL condition and the induction hypothesis, al terms in 
the right hand side of the equations are positive and so (f/j+i)° — (Vj+i)^ > 0. 

□ 

The next Discrete-L^-Contraction property is the analogous of the one 
presented in (Theorem 6.2). 



Theorem 3.6. [L^- Contraction] Under the assumptions of Corollary 3.5 let 
= [(^i)i']^^'" o,nd (yj)i = [(Zj)f]^/''". The next contractions property 

holds, 



I-l 



I-l 



(3.16) 
for all j 



i=i i=i 
., J. 
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Note 2. A mass decay property is a direct consequence of (3.16). 

Proof. We will prove the result for j = 1 for simplicity, but and induction 
method could be used to prove it for a general j. We have the next two 
relations for the solutions at the boundary 



(f/i)° = |^((w^o)j-r(x.)) + /(x.). 



At 
Ax 



{{Z,)]-g^{x;))+g{x,). 



Subtracting both and summing for all i, we get 



i-i n-i i-i \ i-i 

i=l \i=l 1=1 / i=l 

So, if we prove that 

i-i i-i 

i=l i=l 

we are done. To prove this, we first recall that, for any A; = 1, — 1 we 
have 

We call (-ffj)f = {Wj)i—{Zj)^. Note that, thanks to the comparison principle, 
{Hj)^ > 0. Summing up, 

1=1 i=l i=l 1=1 i=l 

I-l 7-1 I~2 I-l 



\k+l 



1=1 i=2 i=l 1=1 



i=l i=l 

Now, remember that iWj)f, {Zj)^, (Hj)^ = thanks to our homogenous 
dirichlet boundary condition. Then, using the above relation with k = K — 1 
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we get 

7-1 7-1 7-1 



(3.17) Y.^H,)f-' = 2j2iH,)r' + iH,)1 + {Hj)t,'-J2iH: 

i=l i=l i=\ 

> E(^.)f-^- 

i=l 

By induction we get then, 

i=l i=l 

for all A; = 0, — 1, which in particular states that, 

1=1 i=l 



K 
3)i 



□ 



As we have said, this contraction property directly implies a total mass 
decay property. 



7-1 7-1 



(3.18) 5^(f/,)° < $^(t/,-i)° 

i=i i=i 

This mass decay is a direct consequence of the homogenous dirichlet bound- 
ary data that we have artificially imposed. We prove now that if we pose the 
scheme in then, the expected conservation of mass property is true. 

Theorem 3.7. [Conservation of mass] Let f G L^(M^) nonnegative. Let 
also (Wj)'^ be the correspondent solution of Problem (3.3) posed in = 
(— oo.oo) X [0, oo) with m > 1. Then 



(3.19) J](t/,)° = 5^(t/,-i)° 

i=i i=i 

for all i = 1, J. 

Proof. Directly from the numerical method we have the next relation for the 
solution in the boundary for every i e Z, 
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(f^i)° = |^[(W^o)J-(lVo)°]+TO°. 



Summing up on alH G Z, 

At 



Ax 



so it is enough to prove that (W^o)i = (^o)° to finish. We can 

i=—oo i=—oo 

subtract from our numerical scheme that for any V/c > 1 we have 
and summing up, 



E (^o)r^ 

i=—oo 

(3.20) 



4 E (^o)' - E (^o)ti - E (^o)ti - E (^o)^' 

i=— oo j=— oo i=— oo i=—oo 

oo oo 

2 E (^o)' - E (^o)^'- 

i=— oo j=— oo 

oo 

Now, lets call = (W^o)i^- Since we have finite total mass, we can 



assume that = 1. Then relation (3.20) can be interpreted as the next 
recurrence succession, 

ak+2 = 2afc+i - ak with ao = 1 

The general solution of this recurrence is = c + kd, where c and d are 
two constants which depends on the initial condition. But here we have only 
the initial condition ao = 1, which gives Ofc = 1 + kd. What we want at this 
point is a second initial conditions ai = 1 and then the general solution will 
be 



E 



(lV„)f = at = l 

Proceeding by contradiction, assume ai < 1, then 

ttk = I + k{ai — 1) — oo 

which is a contraction with the maximum principle. The other options is 
assuming ai > 1 but then, you get ak oo which is again a contradiction 
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with the resuh of loosing mass in the bounded domain. So ai = 1 is the only 
option. Then we get the desired result 

oo oo 

J2 ml = a, = ao=Y. (^o)°- 

i=— oo i=— oo 

□ 



4 Comments 



4.1 Proofs for iV > 1 

As we have said before, all the proofs written here are also valid when N 
is greater than one, but in order to convince the reader we will at least 
formulate the numerical scheme in this case. 

We need to introduce some multi index notation for the spatial discretiza- 
tion, i = (zi, iAf) with ii = 0, // for / = 1, iV where /; is the number 
of nodes of our mesh in the Z-esim dimensional direction. We will also use 
say that 

lz = (0,..., ^,...,0) 

l—esim 



In this way. Problem (3.3) becomes 



(4.1) 

/ N 



1=1 



Ax2 



At 
Ax 



where the solution of 



0, <i < I, 
0<k<K, 

0<i < I, 
otherwise, 



(Wo)', = r (x,), 



0, 



0, <i < I, 
0<k<K, 

<i < I, 

otherwise. 
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is used to start the numerical method. 

With this multi-index notation the proofs for the local truncation error, 
existence, uniqueness and convergence are valid without any change. 



4.2 Convergence to the solution in the whole space 

Until now, we have only proved the convergence of the solution {Wj)^ to the 



numerical method (3.3) to the real solution w of Problem (2.3) posed in the 
bounded domain [X, X] x [0,X]. In the following, we will call (W^)'^ and 
to those solutions. It is proven in [6] that, if we call w to the solution in 
the whole space, then 

II Xii X^co „ 

\\W-W ||l°o(]rJVx[o,T]) — > U. 

But now, if we call E = max \w{xi,yk,tj) — {Wf)^\ then 
E = max \w{xi, yk, tj) - {xi, yk, tj) + (x,,, yk, tj) - {Wf)^] 



< max \w{xi, yk, tj) - w (xj, yk, tj) \ + max \w {xi, yk, tj) - {W^ ) 

i,j,k 



< 



4.3 Required regularity 

In all the previous arguments, we have computed all the needed derivatives 



without taking any care of the regularity of the solution to (2.3). We are 



not going to develop that theory of regularity, but at least we will post the 



required one to make all the proofs valid. In the Theorem 3.1 a taylor 
expansion of the solution is made in several ways. We need the residues 
of the taylor expansions, the higher order derivatives, to exists and to be 
bounded, that is: 

1. w{-, -, t) G C\n) for all t G (0, T). 

2. w{.,0,t) G C^{Td) for all t G (0,T). 

3. w{x, 0, .) G C2((0, T]) for all x G T^. 

The C°° regularity of nonnegative solutions of the FPME with integrable 
initial data is settled in the paper [9], where the case of the whole space is 
treated. For a bounded domain similar results are still to be proved. 
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4.4 Signed initial data 



We have assume that the initial data / is a nonnegative function only for 
simplicity. One can easy observe that all the proofs are valid for / with sign 
with very small changes. Again, problems could come proving the regularity 
required for the theoretical solutions. Lets call 

bmin = min{/"'(x),0}. 



Then, in Theorem 3.2 the same argument holds to prove (Uj)^ > bmin and 
so, the new maximum/minimum principle states, under the same assump- 
tions of the previous one, that for all i,j, k 

bmin — (J-^j^i — ^max- 



In Theorem 3.4 nothing change since we are always talking about errors. 



and so we are working with absolute values. 



5 Avoiding the regularity problems: Com- 
parison with the problem in the whole M^^^ 

As we have said before, in many of our proofs we have assumed enough reg- 
ularity for our theoretical solutions in the bounded domain that aloud us to 
use taylor expansion until the required order. Great part of that regularity 
is proven in [6] but not all. In this section we give another numerical ap- 
proach that avoids this problem but in return gives an extra condition for 
the convergence. 



The idea passes through comparing the numerical solution to (3.3) posted 
in the bounded domain [—X, X] x [0, X] (here appears the first extra condi- 
tion which s ays t hat Y and X have to be proportionals) with the theoretical 



solution to (2.2) posed in M;!^^^. It is proven in [0] that we have the re- 
quired regularity for the solution w. But a new problem comes: w 7^ in F/i 
and so we cannot have convergence for a fixed domain. The idea is making 
f2 — > as Ax — > with a certain rage of convergence. 

In [I2J, an upper bound for the solution with compactly supported initial 



data is found passing through the Barenblatt solutions of problem (2.1 ). The 
upped bound is, 

uUx,t) = it + l)-'^F{\x\it + ir^), 
where F(0 < and 

a = — — -, p 



A^(m + 1) + 1' A^(m + 1) 
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Since —a + fi{N + 1) = /3, we have the next bound in Yy, depending only on 

<,(x, t) < c{t + < c{t + ir^^ <c-{T+ ir^^. 



Then, if we impose the second condition 



(5.1) 



IXI > K 



1 



we can adapt the proofs of Theorems 3A_ and 3^ to obtain the whised con- 
vergence. 



In Theorem 3.1 , the local truncation error in the interior of nodes and 
in Yd still being the same but is not zero anymore in Y^. Now if (xj, yk) E Yh, 



X\ 



and so A = O(Ax^) as before. 

In Theorem |3.4 , again the unique change is that the error in Y^ is not 
zero. But, if {xi,yk) G Yh, 



[W 



Jli 



And so Ej = 0{Ax + At). 
We get then the next result. 



Theorem 5.1. Let w be the solution to Problem (2.2) (posed in M. ) and 
iyVj)^ be the solution to system (3.3) (posed in the bounded domain VL = 
[—X, X] X [0, X]) with m > 1 and compactly supported initial data f G 
L^{R^). Assume that: 

1. There exists two constants C{m, f),D > such that 

DAx <At< C(m, f)Ax. 

2. There exists a constant K > such that 

1 



IXI > K- 



Then 



Ax 



max\w{xi,yk,tj) - {Wj)'l\ = 0{Ax + At). 



Note 3. Condition 2 says that \X 



Ax-5>0 



oo and so Q 
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6 Extension to a more general fractional dif- 
fusion equation 

It is also possible to formulate a numerical method for more general equations 



(6.1) 



du 



X, t) + (-A)^/V(M)(a;, t) = 0, x G M^, t > 



m(x,0) = /(x), 



X G 



a if E C^(]R) such that > and locally bounded. We also need 

existence of >^^^ G C^(]R). In this case, we have an associated extension 
problem 



f Ait;(a;,y,t) = 0, 



(6.2) 



x^W , y > 0, t > 0, 



^^^^(x,0,t) = ^(a;,0,t), a;GM^ y = 0, t>0, 
L w(x,0,0) = <^(/(x)), xgM^, 



as before. In [B] they study the theory for this problem with 

Lp{u) = \0g{u + 1), 

as the natural limit as m — )■ of the problem with 



m 



We can also state the corresponding finite difference method associated to 

this problem posed in the bounded domain, 

(6.3) 



Ax2 



0, Q <i < <k < K 



f At 



= - i^O<^<I 

iWj)i = 0, otherwise, 
where the solution of 



Ax"^ 



(W^o)° = V^(/(x.)), 
I iWo)1 = 0, 



0, 0<i<I,0<k<K 
if < i < J, 
otherwise. 



21 



is used to start the numerical method. 

All proofs can be adapted, without any extra effort, to this method. Of 



course we have to formulate a more general constant (3.6) adapted to (p 



We remember that the restriction for this constant comes from the required 



positivity of the coefficient 1 — V9'(^) At/Ax that appears in (3.9), and so, the 
constant becomes 



(6.4) C{ipJ)=[ max {y,'(x)} 

All the calculus work in the same way using this constant and again the key 
point of convergence pass through note that the function 

x-y 

is Ci((0,oo) X (0 X oo)) . 
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7 Numerical results 



7.1 Analysis of the errors 

We present now two different error analysis of the numerical solutions. The 
first one is obtained comparing solutions with some decreasing Ax with a 
solution generated with a very small Ax. 



1 




50 100 150 200 250 300 350 400 450 500 lo^ io' 

1og(Nodes) 



Figure 1: First analysis of the error for m=l and m=2 
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The second way of computing the errors can only be done for the case 
m = 1 since we have and explicit solution when the initial data is a Dirac 
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delta. The idea is to take as initial data for our method the explicit solution 
given by 

u(x,t) = C^^^^, 

at time t=l. So the error will be be computed with the difference of the 
solution of the method with initial data 



fix) = C 



1 



X P + 1 



at time T = 1 and the real solution at time T = 1, that is, the solution of 
the problem with initial data the dirac Delta at time T = 2. We are going 
to compare the a real solution in the whole space with a numerical solution 
computed in a bounded domain so we choose a very large domain in order 
to minimize the errors that comes from the tails. The chosen domain is 
n = [-100,100] X [0,100]. 




Figure 2: Second error analysis of the error for m=l 
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A third way of computing errors for m 7^ 1 is possible thanks to the 
Barenblatt formula introduced by J.L. Vazquez in [12]. In that paper this 
numerical method is used to compute some Barenblatt profiles as next picture 
shows, 




Figure 3: Computed Barenblatt profiles for m = 1,2, 10 with s = 1/2 
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7.2 Graphics of some solutions 

We next present some graphics with the numerical results obtained with 
initial data 

/(x) =Ce-(i-)W)X[-i,i](a;). 

In Figures 4 and 5 we show the numerical solutions for two small m where 
the expected fat tail is observed. 




Figure 4: Numerical solution for m = 1 




Figure 5: Numerical solution for m = 2 
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In the Figure 6 we observe the typical very slow diffusion of the porous 
medium equation with high values of m. 




-2 



Figure 6: Numerical solution for m = 10 

A case with different diffusion is presented in Figure 7 as an example 
where the numerical method is used for more general ^p. 



0.9-^ 
0.8 




Figure 7: Numerical solution for '^[u) = log{u + 1) 
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